function ll = lifetime_bonacich_log_likelihood(theta,n,in_deg_years,out_deg_years,bonacich_years,startyear,lastyear,censored)

delta = theta(1);
zeta = theta(2);
rho = theta(3:end);

t_max = size(in_deg_years,1);

%% Compute hazard functions.
lifetime_hazard = zeros(n,t_max);
for i=1:n
    for t=1:t_max
        yy=2002+t;
        if startyear(i)<=yy && yy<=lastyear(i)
        %if firm is active at that year, compute a nonzero hazard
            dummy = zeros(length(rho),1);
            for k=1:length(rho)
                dummy(k) = (in_deg_years(t,i)==k-1);
            end
            %lifetime_hazard(i,t) = delta + zeta*1/(1+out_deg_years(t,i)) + rho*dummy;
            lifetime_hazard(i,t) = delta + zeta/bonacich_years(t,i) + rho*dummy;
        end
    end
end
lifetime_hazard = max(lifetime_hazard,0);

%% Compute cumulative hazard functions.
lifetime_cum_hazard = nan(n,t_max);
for i=1:n
    for t=1:t_max
        lifetime_cum_hazard(i,t) = sum(lifetime_hazard(i,1:t));
    end
end

%% Compute log-likelihood.
dummy = 0;
for i=1:n
    t = lastyear(i)-2002;
    dummy = dummy + (1-censored(i))*log(lifetime_hazard(i,t)) - lifetime_cum_hazard(i,t);
end
ll = dummy;
if(isinf(ll))
    ll = -10^10;
end

